In [1]:
import os
import warnings
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import plotly as py
import chart_studio
import plotly.offline
import plotly.graph_objs as go
import cufflinks as cf
from plotly.graph_objs import * 
%matplotlib inline
plotly.offline.init_notebook_mode(connected=True)
cf.go_offline()
warnings.filterwarnings('ignore')
In [2]:
data = pd.read_csv('air_data.csv',encoding='ANSI')
print(data.shape)
data.head()

(62988, 44)
Out[2]:
MEMBER_NO FFP_DATE FIRST_FLIGHT_DATE GENDER FFP_TIER WORK_CITY WORK_PROVINCE WORK_COUNTRY AGE LOAD_TIME ... ADD_Point_SUM Eli_Add_Point_Sum L1Y_ELi_Add_Points Points_Sum L1Y_Points_Sum Ration_L1Y_Flight_Count Ration_P1Y_Flight_Count Ration_P1Y_BPS Ration_L1Y_BPS Point_NotFlight
0 54993 2006/11/2 2008/12/24 男 6 . 北京 CN 31.0 2014/3/31 ... 39992 114452 111100 619760 370211 0.509524 0.490476 0.487221 0.512777 50
1 28065 2007/2/19 2007/8/3 男 6 NaN 北京 CN 42.0 2014/3/31 ... 12000 53288 53288 415768 238410 0.514286 0.485714 0.489289 0.510708 33
2 55106 2007/2/1 2007/8/30 男 6 . 北京 CN 40.0 2014/3/31 ... 15491 55202 51711 406361 233798 0.518519 0.481481 0.481467 0.518530 26
3 21189 2008/8/22 2008/8/23 男 5 Los Angeles CA US 64.0 2014/3/31 ... 0 34890 34890 372204 186100 0.434783 0.565217 0.551722 0.448275 12
4 39546 2009/4/10 2009/4/15 男 6 贵阳 贵州 CN 48.0 2014/3/31 ... 22704 64969 64969 338813 210365 0.532895 0.467105 0.469054 0.530943 39

5 rows × 44 columns

导入数据集csv,以ANSI编码格式打开

一.在接下来的步骤中，将对数据进行前期观测，并尝试着分析某些数据间是否存在着某种关联，同时也对LRFMC模型中需要用到的属性，进行观测。

In [3]:
plt.subplots(figsize=(14,14))
sns.heatmap(data.corr(),annot=True)
Out[3]:
<AxesSubplot:>

使用heatmap大体查看数据各列间是否存在显著联系，没看出什么

In [4]:
sns.pairplot(data=data)
Out[4]:
<seaborn.axisgrid.PairGrid at 0x1c9ce562f50>

使用pairplot查看数据分布之间的关联

In [5]:
data.GENDER.value_counts()
Out[5]:
男    48134
女    14851
Name: GENDER, dtype: int64

查看乘客性别分布数

In [6]:
data.WORK_PROVINCE.value_counts()
Out[6]:
广东               17507
北京                8014
上海                4994
辽宁                4182
新疆                2512
                 ...  
廣西省                  1
AICHIKEN             1
SOUTU HOLLAND        1
新潟县                  1
JAWATIMUR            1
Name: WORK_PROVINCE, Length: 1183, dtype: int64

查看乘客工作省份分布数

In [7]:
columns = [x for x in data.columns]
percent=[]
for kolom in columns:
    percent.append(round(data[kolom].isnull().sum()/data[kolom].shape[0]*100, 2))
    
explore = data.describe(percentiles = [], include = 'all').T 
explore['missing'] = len(data) - explore['count'] 
explore['%'] = percent
explore = explore[['missing','%','min','max']]
explore = explore.replace(np.nan, '-', regex=True)
data = data[data['SUM_YR_1'].notnull()]
data = data[data['SUM_YR_2'].notnull()]
print(explore)
data

                         missing     %    min        max
MEMBER_NO                    0.0  0.00    1.0    62988.0
FFP_DATE                     0.0  0.00      -          -
FIRST_FLIGHT_DATE            0.0  0.00      -          -
GENDER                       3.0  0.00      -          -
FFP_TIER                     0.0  0.00    4.0        6.0
WORK_CITY                 2269.0  3.60      -          -
WORK_PROVINCE             3248.0  5.16      -          -
WORK_COUNTRY                26.0  0.04      -          -
AGE                        420.0  0.67    6.0      110.0
LOAD_TIME                    0.0  0.00      -          -
FLIGHT_COUNT                 0.0  0.00    2.0      213.0
BP_SUM                       0.0  0.00    0.0   505308.0
EP_SUM_YR_1                  0.0  0.00    0.0        0.0
EP_SUM_YR_2                  0.0  0.00    0.0    74460.0
SUM_YR_1                   551.0  0.87    0.0   239560.0
SUM_YR_2                   138.0  0.22    0.0   234188.0
SEG_KM_SUM                   0.0  0.00  368.0   580717.0
WEIGHTED_SEG_KM              0.0  0.00    0.0  558440.14
LAST_FLIGHT_DATE             0.0  0.00      -          -
AVG_FLIGHT_COUNT             0.0  0.00   0.25     26.625
AVG_BP_SUM                   0.0  0.00    0.0    63163.5
BEGIN_TO_FIRST               0.0  0.00    0.0      729.0
LAST_TO_END                  0.0  0.00    1.0      731.0
AVG_INTERVAL                 0.0  0.00    0.0      728.0
MAX_INTERVAL                 0.0  0.00    0.0      728.0
ADD_POINTS_SUM_YR_1          0.0  0.00    0.0   600000.0
ADD_POINTS_SUM_YR_2          0.0  0.00    0.0   728282.0
EXCHANGE_COUNT               0.0  0.00    0.0       46.0
avg_discount                 0.0  0.00    0.0        1.5
P1Y_Flight_Count             0.0  0.00    0.0      118.0
L1Y_Flight_Count             0.0  0.00    0.0      111.0
P1Y_BP_SUM                   0.0  0.00    0.0   246197.0
L1Y_BP_SUM                   0.0  0.00    0.0   259111.0
EP_SUM                       0.0  0.00    0.0    74460.0
ADD_Point_SUM                0.0  0.00    0.0   984938.0
Eli_Add_Point_Sum            0.0  0.00    0.0   984938.0
L1Y_ELi_Add_Points           0.0  0.00    0.0   728282.0
Points_Sum                   0.0  0.00    0.0   985572.0
L1Y_Points_Sum               0.0  0.00    0.0   728282.0
Ration_L1Y_Flight_Count      0.0  0.00    0.0        1.0
Ration_P1Y_Flight_Count      0.0  0.00    0.0        1.0
Ration_P1Y_BPS               0.0  0.00    0.0   0.999989
Ration_L1Y_BPS               0.0  0.00    0.0   0.999993
Point_NotFlight              0.0  0.00    0.0      140.0
Out[7]:
MEMBER_NO FFP_DATE FIRST_FLIGHT_DATE GENDER FFP_TIER WORK_CITY WORK_PROVINCE WORK_COUNTRY AGE LOAD_TIME ... ADD_Point_SUM Eli_Add_Point_Sum L1Y_ELi_Add_Points Points_Sum L1Y_Points_Sum Ration_L1Y_Flight_Count Ration_P1Y_Flight_Count Ration_P1Y_BPS Ration_L1Y_BPS Point_NotFlight
0 54993 2006/11/2 2008/12/24 男 6 . 北京 CN 31.0 2014/3/31 ... 39992 114452 111100 619760 370211 0.509524 0.490476 0.487221 0.512777 50
1 28065 2007/2/19 2007/8/3 男 6 NaN 北京 CN 42.0 2014/3/31 ... 12000 53288 53288 415768 238410 0.514286 0.485714 0.489289 0.510708 33
2 55106 2007/2/1 2007/8/30 男 6 . 北京 CN 40.0 2014/3/31 ... 15491 55202 51711 406361 233798 0.518519 0.481481 0.481467 0.518530 26
3 21189 2008/8/22 2008/8/23 男 5 Los Angeles CA US 64.0 2014/3/31 ... 0 34890 34890 372204 186100 0.434783 0.565217 0.551722 0.448275 12
4 39546 2009/4/10 2009/4/15 男 6 贵阳 贵州 CN 48.0 2014/3/31 ... 22704 64969 64969 338813 210365 0.532895 0.467105 0.469054 0.530943 39
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
62982 16415 2013/1/20 2013/1/20 女 4 北京 . CN 35.0 2014/3/31 ... 0 0 0 0 0 0.000000 1.000000 0.000000 0.000000 0
62983 18375 2011/5/20 2013/6/5 女 4 广州 广东 CN 25.0 2014/3/31 ... 12318 12318 12123 12318 12123 1.000000 0.000000 0.000000 0.000000 22
62984 36041 2010/3/8 2013/9/14 男 4 佛山 广东 CN 38.0 2014/3/31 ... 106972 106972 56506 106972 56506 1.000000 0.000000 0.000000 0.000000 43
62985 45690 2006/3/30 2006/12/2 女 4 广州 广东 CN 43.0 2014/3/31 ... 0 0 0 0 0 1.000000 0.000000 0.000000 0.000000 0
62986 61027 2013/2/6 2013/2/14 女 4 广州 广东 CN 36.0 2014/3/31 ... 0 0 0 0 0 0.000000 1.000000 0.000000 0.000000 0

62299 rows × 44 columns

统计数据缺失数量、确实比例、确实最大最小值

In [8]:
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['AGE'], ax=ax_box)
sns.distplot(data['AGE'], ax=ax_hist)
plt.show()

查看AGE分布情况

In [9]:
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['LAST_TO_END'], ax=ax_box)
sns.distplot(data['LAST_TO_END'], ax=ax_hist)
plt.show()

查看LAST_TO_END分布情况

In [10]:
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['FLIGHT_COUNT'], ax=ax_box)
sns.distplot(data['FLIGHT_COUNT'], ax=ax_hist)
plt.show()

查看FLIGHT_COUNT分布情况

In [11]:
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['SEG_KM_SUM'], ax=ax_box)
sns.distplot(data['SEG_KM_SUM'], ax=ax_hist)
plt.show()

查看SEG_KM_SUM分布情况

In [12]:
data['FFP_DATE']=pd.to_datetime(data['FFP_DATE'])
data['LOAD_TIME']=pd.to_datetime(data['LOAD_TIME'])
from datetime import datetime
def interval_time(dd):
    return dd.days
data['LENGTH']=np.abs(data['LOAD_TIME']-data['FFP_DATE'])
data['LENGTH']=data['LENGTH'].apply(interval_time)
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['LENGTH'], ax=ax_box)
sns.distplot(data['LENGTH'], ax=ax_hist)
plt.show()

将FFP_DATE和LOAD_TIME使用pandas转换为datatime格式，建立新的列LENGTH为LOAD_TIME-FFP_DATE，并查看LENGTH分布情况

In [13]:
f, (ax_box, ax_hist) = plt.subplots(2, sharex=True, gridspec_kw= {"height_ratios": (0.2, 1)},figsize=(10,6))
sns.boxplot(data['avg_discount'], ax=ax_box)
sns.distplot(data['avg_discount'], ax=ax_hist)
plt.show()

查看avg_discount分布情况

In [14]:
data['year'] = data.FFP_DATE.dt.year
data['month'] = data.FFP_DATE.dt.month
data['months'] = data.year.astype(str) + '-' + data.month.astype(str)
data['time']=1
ads = data.groupby('months').time.sum().sort_index()
ads.index = pd.to_datetime(ads.index)
ads.sort_index().ta_plot(
    title='ads', 
    study='sma', periods=[3, 6], study_colors=['lightblue', 'blue'], 
    interpolation='hv', theme='solar'
    )

虽然本数据集中数据本身就只有year和month属性了，不过还是根据我的习惯将FFP_DATE中的year和month提取出来，并强制转换为str格式然后连接，作为下面要用的months。然后在data里添加新列time，作为统计某种数据出现次数的计算量。然后我们将按照months对数据进行分组并计总次数。再将此数据进行索引排列。然后同时使用sma进行以每3和每6月的数据预测。

In [15]:
data['sumc']=data['SUM_YR_1']+data['SUM_YR_2']
adc = data.groupby('months').sumc.sum().sort_index()
adc.index = pd.to_datetime(adc.index)
adc.sort_index().ta_plot(
    title='adc',
    study='sma', periods=[3, 6], study_colors=['lightblue', 'blue'], 
    interpolation='hv', theme='solar',
    vspan={'x0':'2009-02', 'x1':'2010-02', 'color':'lightblue', 'fill':True, 'opacity':.1}
    )

对于数据集中存在的销售数据SUM_YR_1和SUM_YR_2，我们进行求和，并输入到新建的列sumc中，作为后续使用。然后我们将按照months对数据进行分组并求销售和。再将此数据进行索引排列。然后同时使用sma进行以每3和每6月的数据预测。并且在上面的计次中，我们观测到，在2009年的新增用户，大大减少，所以我们猜测当年销售额也将减少，故我们使用vspan框出09年2月到10年2月的数据着重显示。

In [16]:
cts = data.WORK_PROVINCE.value_counts().head(20).reset_index().rename({'index': 'WORK_PROVINCE', 'WORK_PROVINCE': 'Counts'}, axis=1)
argpricebe = data.groupby('WORK_PROVINCE').sumc.sum() / data.WORK_PROVINCE.value_counts()
argprice = argpricebe.head(20)
fig = py.subplots.make_subplots(1, 2, subplot_titles=['cts', 'argprice'])
fig.append_trace(go.Bar(x=cts.WORK_PROVINCE, y=cts.Counts, name='cts'), 1, 1)
fig.append_trace(go.Bar(x=argprice.index, y=argprice.values, name='prs'), 1, 2)
fig.update_layout({'template': 'plotly_dark'})
fig.iplot(legend=True)

在此，我们首先对工作省份进行计次统计，并使用最多的20个省份，同时重置索引。然后对工作省份分类，并将每个省份的sumc除以该省份的出现次数，此举旨在获得该省的用户'单价'。然后我们对于工作省份以及'单价'最多最高的20个省进行可视化。发现由于原数据原因，有一些地名是乱码，不过我们可以不做可视化，只将分类后数据输出，提交给人工分析使用。但是要注意，这里的'单价'，由于数据是每个用户的个人数据为一行，那么此处'单价'的意义为该区域，用户的平均消费水平，而不是用户每次坐飞机的单价。

In [17]:
def unpack_months(field, func, first_group='WORK_PROVINCE'):
    group = data.groupby([first_group, 'months'])
    names = list(data[first_group].value_counts().index)
    tmp = group[field].apply(func).astype(int)
    ret = pd.DataFrame(columns=names)
    for name in names:
        ret[name] = tmp[name]
    ret.index = pd.to_datetime(ret.index)
    return ret.sort_index()

在这里我使用了分组聚合方式，旨在通过日期序列，来进行某特定两项数据的可视化，此unpack_months方式也是我在可视化课程设计中所使用的，受益匪浅。 首先根据原有参数WORK_PROVINCE或后期传入参数，与months进行groupby。然后对于first_group变量进行计次统计并按索引输出为名字列表。 注意tmp中这里的group不是聚合方式，而是上文所设置的df存储，这里使用group中的field数据列，对field使用func，并将结果强行转换格式为int。 而后设置新的dataframe名为ret，并将之前所提取的names输入进ret的columns。 最后根据names中存储的名字，将ret中所存储的数据与tmp所存储的数据进行一一对应。 然后我们设置ret的索引为按照日期，然后返回他可用于外部输出。

In [18]:
adsm = unpack_months('time', sum)
adsm = adsm.iloc[:,:20]
adsm.iplot(
    title='adsm',
    xTitle='Months', yTitle='ads', 
    fill=True, theme='solar', interpolation='hv'
    )

这里我们使用了上方所定义的分组聚合方式，我们输入time和sum，也就是根据工作省份，统计出现time求和的次数，并按照日期展示。当然我们为了展示方便，值提取了前20个最多的省来查看。

In [19]:
adcm = unpack_months('sumc', sum)
adcm = adcm.iloc[:,:20]
adcm.iplot(
    title='adcm',
    xTitle='Months', yTitle='adcm', 
    fill=True, theme='solar', interpolation='hv'
    )

这里我们依然是分组聚合，对于不同城市，按照日期展示sumc。

In [20]:
countries = data.groupby(['WORK_CITY', 'WORK_PROVINCE'])['time', 'sumc'].sum().reset_index()
countries['argpricebe'] = pd.Series(countries.sumc / countries.time)
top_sales = countries.sort_values('time', ascending=False).head(20).sort_index()
colormap = dict(zip(top_sales.WORK_CITY.value_counts().index, py.tools.DEFAULT_PLOTLY_COLORS))
top_sales['colors'] = top_sales.WORK_CITY.map(colormap)
top_sales.iplot(
    kind='bubble', 
    x='WORK_PROVINCE', y='argpricebe', size='time', 
    title='city-prov-sumc-argp',
    yTitle='argpricebe',
    theme='solar', colors=list(top_sales.colors)
    )

然后我们想要查看，对于不同省份的不同城市，他的出现次数，以及用户'单价'。这里的y轴属性表示用户'单价'，圆圈的大小表示出现次数。当然，由于原始数据的地名有些杂乱，我们还无法精准查看。而这也并不是我们课程设计的主要目的，所以我们依然是可以将这份数据输出到csv供专家使用，或着使用tableau快速显示。

二.在上述的步骤中，我们已经完成了能想到的一些前期可视化分析，接下来我们进入到使用k-means针对特定属性，聚类客户。

In [21]:
plt.subplots(figsize=(12,12))
col_n = ['AGE','LENGTH','LAST_TO_END','SEG_KM_SUM','avg_discount','sumc']
data1 = pd.DataFrame(data,columns = col_n)
sns.heatmap(data1.corr(),annot=True)
Out[21]:
<AxesSubplot:>

在这里我们根据LRFMC模型所需要的几个属性，以及一个我个人比较感兴趣的属性AGE进行heatmap，当然也可以添加性别进去，这点只需要将'男'和'女'转变成例如'1'和'0'这样的可输入数据即可。我认为性别和年龄，对于不同聚类类型的用户，也是有其影响的。当然我们后续的工作并不包含这两种属性，但是可以添加进去尝试。

In [22]:
index1=data['SUM_YR_1'] !=0
index2=data['SUM_YR_2'] !=0
index3=(data['SEG_KM_SUM']> 0) & (data['avg_discount'] != 0)  
data=data[(index1 | index2) & index3]
df_lrfmc=data[['MEMBER_NO','LENGTH','LAST_TO_END','FLIGHT_COUNT','SEG_KM_SUM','avg_discount']]
df_lrfmc.columns = ['MEMBER_NO','L','R','F','M','C']
df_lrfmc_drop=df_lrfmc[['L','R','F','M','C']]
df_lrfmc_drop.head()
Out[22]:
L R F M C
0 2706 1 210 580717 0.961639
1 2597 7 140 293678 1.252314
2 2615 11 135 283712 1.254676
3 2047 97 23 281336 1.090870
4 1816 5 152 309928 0.970658

这里我们为了建模的精准性，只保留票价非零的，或者平均折扣率不为0且总飞行公里数大于0的记录。然后将这些数据只提取出我们需要的LRFMC以及用户ID作为df_lrfmc留存，然后我们仅提取出LRFMC保存在df_lrfmc_drop。

In [23]:
from sklearn.preprocessing import StandardScaler
scale = StandardScaler()
dat=df_lrfmc_drop.copy()
kolom_all = [x for x in dat.columns]
for kolom in kolom_all:
      dat[kolom] = scale.fit_transform(np.array(dat[kolom]).reshape(-1,1))
dat
Out[23]:
L R F M C
0 1.435719 -0.944955 14.034129 26.761370 1.295551
1 1.307162 -0.911902 9.073286 13.126970 2.868199
2 1.328392 -0.889866 8.718940 12.653583 2.880973
3 0.658481 -0.416102 0.781591 12.540723 1.994730
4 0.386035 -0.922920 9.923716 13.898848 1.344346
... ... ... ... ... ...
62974 2.076144 -0.460173 -0.706662 -0.805303 -0.065898
62975 0.557051 -0.283888 -0.706662 -0.805303 -0.282311
62976 -0.149422 -0.735617 -0.706662 -0.772338 -2.689906
62977 -1.206183 1.605662 -0.706662 -0.779843 -2.554648
62978 -0.479660 0.603044 -0.706662 -0.786683 -2.392338

62044 rows × 5 columns

使用StandardScaler对数据进行标准化

In [24]:
from sklearn.cluster import KMeans
ks = range(1,11)
inertias=[]
for k in ks :
    kc = KMeans(n_clusters=k,random_state=142)
    kc.fit(dat)
    cluster = kc.fit_predict(dat)
    inertias.append(kc.inertia_)
f, ax = plt.subplots(figsize=(10, 6))
plt.plot(ks, inertias, '-o')
plt.xlabel('Number of clusters, k')
plt.ylabel('Inertia')
plt.xticks(ks)
plt.style.use('ggplot')
plt.title('Choose the Best Number for KMeans')
plt.show()

按照步骤进行k-means的情况，我们先假定不预先知道要分为几类，所以我们要先判断聚类中心数。首先根据手肘法，观察是否在哪个点，数据明显变得平稳，发现在3-6的范围里。

In [25]:
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
for n_clusters in range(3,7):
    clusterer = KMeans(n_clusters=n_clusters,init='k-means++',max_iter=300, random_state=142).fit(dat)
    preds = clusterer.predict(dat)
    score = silhouette_score(dat, preds, metric='mahalanobis')
    print ("For n_clusters = {}. The average silhouette_score is : {}".format(n_clusters, score))

For n_clusters = 3. The average silhouette_score is : 0.23336219118459178
For n_clusters = 4. The average silhouette_score is : 0.2394239225670819
For n_clusters = 5. The average silhouette_score is : 0.24702937931545926
For n_clusters = 6. The average silhouette_score is : 0.24387707446270057

然后我们根据silhouette_score在3到6的范围里评估，最终发现取k为5更合理一点。

In [26]:
from sklearn.cluster import KMeans
kc = KMeans(init='k-means++',n_clusters= 5,random_state=142)
kc.fit(dat)
cluster_labels = kc.labels_
data_c = dat.assign(K_Cluster = cluster_labels)
data_label=dat.assign(K_Cluster = cluster_labels)
data_f=data_c.groupby('K_Cluster').mean()
data_f['count']=pd.value_counts(data_c['K_Cluster'])
data_f.to_csv("data_f.csv")
data_f
Out[26]:
L R F M C count
K_Cluster
0 0.483423 -0.799418 2.482730 2.424282 0.308393 5338
1 -0.700203 -0.414931 -0.160971 -0.160819 -0.255135 24648
2 1.160875 -0.377025 -0.087289 -0.095095 -0.156745 15734
3 0.051942 -0.002219 -0.227343 -0.232068 2.189144 4206
4 -0.314044 1.686414 -0.573989 -0.536776 -0.173210 12118

根据k-means算法，对数据进行聚类，并统计每类的个数。

In [27]:
LRFM=['L','R','F','M','C']
def dist_list(lst):
    plt.figure(figsize=[len(lst)*5,2],)
    i = 1
    cl = ['bgrcmykw'[c] for c in range(len('bgrcmykw'))]
    for col in lst:
        ax = plt.subplot(1,len(lst),i)
        g = data_c.groupby('K_Cluster')
        x = g[col].mean().index
        y = g[col].mean().values
        ax.barh(x,y,color=cl[i-1])
        plt.title(col)
        i = i+1
        
dist_list(LRFM)

('LENGTH','LAST_TO_END','FLIGHT_COUNT','SEG_KM_SUM','avg_discount') 观测5个属性下，优势及劣势的客户群。

In [28]:
x=[1,2,3,4,5]
colors=['blue','green','red','blue','black']
for i in range(5):
    plt.plot(x,kc.cluster_centers_[i],label=('customer%d'%(i)),linewidth=2,color=colors[i],marker='o')
    plt.legend()
plt.xlabel('L R F M C')
plt.ylabel('values')
plt.show()

观测5种客户群，对于5个属性的优势及劣势。

In [29]:
data
Out[29]:
MEMBER_NO FFP_DATE FIRST_FLIGHT_DATE GENDER FFP_TIER WORK_CITY WORK_PROVINCE WORK_COUNTRY AGE LOAD_TIME ... Ration_P1Y_Flight_Count Ration_P1Y_BPS Ration_L1Y_BPS Point_NotFlight LENGTH year month months time sumc
0 54993 2006-11-02 2008/12/24 男 6 . 北京 CN 31.0 2014-03-31 ... 0.490476 0.487221 0.512777 50 2706 2006 11 2006-11 1 473748.0
1 28065 2007-02-19 2007/8/3 男 6 NaN 北京 CN 42.0 2014-03-31 ... 0.485714 0.489289 0.510708 33 2597 2007 2 2007-2 1 338917.0
2 55106 2007-02-01 2007/8/30 男 6 . 北京 CN 40.0 2014-03-31 ... 0.481481 0.481467 0.518530 26 2615 2007 2 2007-2 1 328600.0
3 21189 2008-08-22 2008/8/23 男 5 Los Angeles CA US 64.0 2014-03-31 ... 0.565217 0.551722 0.448275 12 2047 2008 8 2008-8 1 241850.0
4 39546 2009-04-10 2009/4/15 男 6 贵阳 贵州 CN 48.0 2014-03-31 ... 0.467105 0.469054 0.530943 39 1816 2009 4 2009-4 1 255262.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
62974 11163 2005-05-08 2005/8/26 男 4 NaN NaN CN 34.0 2014-03-31 ... 0.000000 0.000000 0.997506 1 3249 2005 5 2005-5 1 960.0
62975 30765 2008-11-16 2013/11/30 男 4 TAIPEI NaN TW 38.0 2014-03-31 ... 0.000000 0.000000 0.997506 0 1961 2008 11 2008-11 1 910.0
62976 10380 2010-07-08 2011/6/21 男 4 贵阳市 贵州省 CN 33.0 2014-03-31 ... 0.500000 0.995327 0.000000 1 1362 2010 7 2010-7 1 284.0
62977 16372 2012-12-20 2012/12/20 男 4 桃园 NaN TW 47.0 2014-03-31 ... 1.000000 0.000000 0.000000 0 466 2012 12 2012-12 1 330.0
62978 22761 2011-04-14 2011/4/14 男 4 汕头 广东省 CN 48.0 2014-03-31 ... 0.000000 0.000000 0.000000 0 1082 2011 4 2011-4 1 370.0

62044 rows × 50 columns

In [30]:
data['K_Cluster']=cluster_labels
sns.scatterplot(x=data['AGE'],y=data['FLIGHT_COUNT'],hue=data['K_Cluster'])
Out[30]:
<AxesSubplot:xlabel='AGE', ylabel='FLIGHT_COUNT'>

根据客户群体，查看年龄与航班次数的关系

In [31]:
sns.scatterplot(x=data['AGE'],y=data['sumc'],hue=data['K_Cluster'])
Out[31]:
<AxesSubplot:xlabel='AGE', ylabel='sumc'>

根据客户群体，查看年龄与消费的关系，这证明了我们之前猜测的，年龄确实在某种程度上，与不同的消费群体有某种关系。相信性别也有其影响，可以像前文我所描述的那样，将'男'和'女'转变成例如'1'和'0'这样的可输入数据，然后就可以进行类似的观测，当然我们也可以将这两样数据放入模型训练中，以获得更精准的群体关系。

In [32]:
ba = pd.read_csv("data_f.csv",encoding='utf-8',engine='python')
ba.iplot(kind = 'bar',x = 'K_Cluster',y = 'count')

不同类别用户群体的数量情况

In [33]:
data_c = data_c[~data_c.index.duplicated()]
sns.pairplot(data=data_c, hue='K_Cluster')
Out[33]:
<seaborn.axisgrid.PairGrid at 0x1caab7547c0>

上图展示了不同客户群体，对于LRFMC五种属性的分布关系。例如我们可以看出，对于客户群体0，在FM关系，或FC关系，则有明显的关联。而对于那些杂乱的图表，我认为他们确实是无关系的杂乱的，不过我们也可以分别提取5种客户群体，依次查看分布。

In [34]:
categories = ['L','R','F','M','C']
fig = go.Figure()
fig.add_trace(go.Scatterpolar(
      r=ba.loc[0,"L":"C"],
      theta=categories,
      fill='toself',
      name='A'
))
fig.add_trace(go.Scatterpolar(
      r=ba.loc[1,"L":"C"],
      theta=categories,
      fill='toself',
      name='B'
))
fig.add_trace(go.Scatterpolar(
      r=ba.loc[2,"L":"C"],
      theta=categories,
      fill='toself',
      name='C'
))
fig.add_trace(go.Scatterpolar(
      r=ba.loc[3,"L":"C"],
      theta=categories,
      fill='toself',
      name='D'
))
fig.add_trace(go.Scatterpolar(
      r=ba.loc[4,"L":"C"],
      theta=categories,
      fill='toself',
      name='E'
))
fig.update_layout(
  polar=dict(
    radialaxis=dict(
      visible=True,
      range=[-2, 3]
    )),
  showlegend=True
)

fig.show()

这可以明显的看出5种客户群体，在LRFMC不同属性上的行为分布

In [35]:
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
n_clusters = 5
sklearn_pca = PCA(n_components = 2)
Y_sklearn = sklearn_pca.fit_transform(dat)
kmeans = KMeans(n_clusters= n_clusters, max_iter=600,random_state=150, algorithm = 'auto')
%time fitted = kmeans.fit(Y_sklearn)
prediction = kmeans.predict(Y_sklearn)
df_pred=data.assign(Pre_K_Cluster = prediction)
targets = [0,1,2,3,4]
plt.figure(figsize = (10,8))
plt.scatter(Y_sklearn[:, 0], Y_sklearn[:, 1],c=prediction ,s=50, cmap='viridis')
centers2 = fitted.cluster_centers_
plt.scatter(centers2[:, 0], centers2[:, 1],c='red', s=300, alpha=0.6);

CPU times: total: 6.98 s
Wall time: 676 ms

这是我的另一种思考和尝试，即利用PCA将数据降至2维到Y_sklearn保存，再利用kmeans的pridict方法来将这些Y_sklearn进行分类，并返回标签到原数据上。

In [36]:
df_pred
Out[36]:
MEMBER_NO FFP_DATE FIRST_FLIGHT_DATE GENDER FFP_TIER WORK_CITY WORK_PROVINCE WORK_COUNTRY AGE LOAD_TIME ... Ration_L1Y_BPS Point_NotFlight LENGTH year month months time sumc K_Cluster Pre_K_Cluster
0 54993 2006-11-02 2008/12/24 男 6 . 北京 CN 31.0 2014-03-31 ... 0.512777 50 2706 2006 11 2006-11 1 473748.0 0 3
1 28065 2007-02-19 2007/8/3 男 6 NaN 北京 CN 42.0 2014-03-31 ... 0.510708 33 2597 2007 2 2007-2 1 338917.0 0 3
2 55106 2007-02-01 2007/8/30 男 6 . 北京 CN 40.0 2014-03-31 ... 0.518530 26 2615 2007 2 2007-2 1 328600.0 0 3
3 21189 2008-08-22 2008/8/23 男 5 Los Angeles CA US 64.0 2014-03-31 ... 0.448275 12 2047 2008 8 2008-8 1 241850.0 0 3
4 39546 2009-04-10 2009/4/15 男 6 贵阳 贵州 CN 48.0 2014-03-31 ... 0.530943 39 1816 2009 4 2009-4 1 255262.0 0 3
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
62974 11163 2005-05-08 2005/8/26 男 4 NaN NaN CN 34.0 2014-03-31 ... 0.997506 1 3249 2005 5 2005-5 1 960.0 2 1
62975 30765 2008-11-16 2013/11/30 男 4 TAIPEI NaN TW 38.0 2014-03-31 ... 0.997506 0 1961 2008 11 2008-11 1 910.0 2 1
62976 10380 2010-07-08 2011/6/21 男 4 贵阳市 贵州省 CN 33.0 2014-03-31 ... 0.000000 1 1362 2010 7 2010-7 1 284.0 1 4
62977 16372 2012-12-20 2012/12/20 男 4 桃园 NaN TW 47.0 2014-03-31 ... 0.000000 0 466 2012 12 2012-12 1 330.0 4 4
62978 22761 2011-04-14 2011/4/14 男 4 汕头 广东省 CN 48.0 2014-03-31 ... 0.000000 0 1082 2011 4 2011-4 1 370.0 4 4

62044 rows × 52 columns

In [37]:
data_c
Out[37]:
L R F M C K_Cluster
0 1.435719 -0.944955 14.034129 26.761370 1.295551 0
1 1.307162 -0.911902 9.073286 13.126970 2.868199 0
2 1.328392 -0.889866 8.718940 12.653583 2.880973 0
3 0.658481 -0.416102 0.781591 12.540723 1.994730 0
4 0.386035 -0.922920 9.923716 13.898848 1.344346 0
... ... ... ... ... ... ...
62974 2.076144 -0.460173 -0.706662 -0.805303 -0.065898 2
62975 0.557051 -0.283888 -0.706662 -0.805303 -0.282311 2
62976 -0.149422 -0.735617 -0.706662 -0.772338 -2.689906 1
62977 -1.206183 1.605662 -0.706662 -0.779843 -2.554648 4
62978 -0.479660 0.603044 -0.706662 -0.786683 -2.392338 4

62044 rows × 6 columns

In [38]:
df_pred['K_Cluster'].value_counts()
Out[38]:
1    24648
2    15734
4    12118
0     5338
3     4206
Name: K_Cluster, dtype: int64
In [39]:
df_pred['Pre_K_Cluster'].value_counts()
Out[39]:
4    20295
1    18631
2    15933
0     3913
3     3272
Name: Pre_K_Cluster, dtype: int64

我们比较了两种方式的最终聚类结果，发现是大同小异，两者的误差我肉眼观测大约在15%？

在这里我们的项目结束了，我仍在思考是否能使用其他来解决类似这样格式的数据集，例如GAN，CNN，又或是使用one hot encoding，pytorch等来帮助我，目前我已经有了一些头绪关于使用pytorch来完成它。如果后续我搞懂了如何实现，也会添加到这后面的。